#include "../include/Matriz.h"
using namespace std;

Matriz::Matriz(int n, int m){
    matriz.clear();
    alto = n;
    ancho = m;
    vector<TFloat> p;
    precision = 52;
    TFloat t(precision);
    t = 0;
    for (int i = 0; i < ancho; i++)
        p.push_back(t);
    for (int i = 0; i < alto; i++)
        matriz.push_back(p);
}

Matriz::Matriz (int n){
    TFloat t(52);
    t=1;
    Matriz a(n,n);
    for (int i=0; i<n; i++){
        for (int j=0; j<n; j++){
            if (i==j) a[i][j] = t;
        }
    }
    alto = n;
    ancho = n;
    matriz = a.matriz;
}

Matriz::Matriz (int n, int m, int prec){
    matriz.clear();
    alto = n;
    ancho = m;
    vector<TFloat> p;
    precision = (size_t) prec;
    TFloat t(precision);
    t = 0;
    for (int i = 0; i < ancho; i++)
        p.push_back(t);
    for (int i = 0; i < alto; i++)
        matriz.push_back(p);
}

Matriz::Matriz (vector<TFloat> v){
    Matriz res(v.size(),1);
    for (uint i=0; i<v.size(); i++){
        res[i][0] = v[i];
    }
    alto = res.alto;
    ancho = res.ancho;
    matriz = res.matriz;
}

void Matriz::cargar(std::istream &is, int pre){
    precision = (size_t) pre;
	matriz.clear();
	vector<TFloat> p;
	is >> alto;
	is >> ancho;
	for (int i = 0; i < alto; i++){
		for (int j = 0; j < ancho; j++){
		    TFloat t (precision);
			is >> t;
			p.push_back(t);
		}
		matriz.push_back(p);
		p.clear();
	}
}

void Matriz::mostrar(std::ostream &os){
    os.precision(5);
    os.setf(ios::fixed,ios::floatfield);
	for (int i=0 ; i<alto; i++){
		os << "[";
		for (int j=0; j<ancho; j++){
		    if (matriz[i][j].dbl() < 0) os.precision(4);
		    else os.precision(5);
		    if (j == ancho-1) os << matriz[i][j];
			else os << matriz[i][j] << ",";
		}
        os << "]" << endl;
	}
}

size_t min(size_t a, size_t b){
    if (a < b)
        return a;
    else
        return b;
}

Matriz Matriz::multiplicar(Matriz a){
    if (ancho!=(a.alto)) cerr << "Estás tratando de multiplicar matrices incompatibles" << endl;
    Matriz b (alto, (a.ancho), min(precision, a.precision));
    for (int i = 0; i < alto; i++){
        for(int j = 0; j < a.ancho; j++){
            b[i][j] = 0;
            for(int k = 0; k < ancho; k++){
                b[i][j] = b[i][j] + matriz[i][k].dbl() * a[k][j].dbl();
            }
        }
    }
    return b;
}

Matriz Matriz::multiplicar(TFloat t){
    Matriz res(alto, ancho, precision);
    for (int i = 0; i < alto; i++)
        for (int j = 0; j < ancho; j++)
            res[i][j] = matriz[i][j]*t;
    return res;
}

Matriz Matriz::sumar(Matriz otra){
    if (alto != otra.alto or ancho != otra.ancho) cerr << "Estás tratando de sumar matrices incompatibles" << endl;
    Matriz res(alto, ancho, min(precision, otra.precision));
    for (int i = 0; i < alto; i++)
        for (int j = 0; j < ancho; j++)
            res[i][j] = matriz[i][j] + otra[i][j];
    return res;
}

Matriz Matriz::restar(Matriz otra){
    if (alto != otra.alto or ancho != otra.ancho) cerr << "Estás tratando de restar matrices incompatibles" << endl;
    Matriz res(alto, ancho, min(precision, otra.precision));
    for (int i = 0; i < alto; i++)
        for (int j = 0; j < ancho; j++)
            res[i][j] = matriz[i][j] - otra[i][j];
    return res;
}

Matriz Matriz::trasponer(){
    Matriz b (ancho, alto, precision);
    for (int i=0; i<alto; i++){
        for(int j=0; j<ancho; j++){
            b[j][i] = matriz[i][j];
        }
    }
    return b;
}

vector<TFloat> Matriz::diagonal(){
    vector<TFloat> diag;
    for(int i=0; i<alto; i++){
        for (int j=0; j<ancho; j++){
            if (i==j) diag.push_back(matriz[i][j]);
        }
    }
    return diag;
}

vector<TFloat> restar(vector<TFloat> a, vector<TFloat> b){
    vector<TFloat> res (a.size());
    if (a.size()!=b.size()) cerr<<"Estas tratando de restar dos vectores incompatibles";
    for (uint i=0; i<a.size(); i++){
        res[i] = a[i] - b[i];
    }
    return res;
}

Matriz Matriz::diagonal(vector<TFloat> v){
    Matriz res(v.size(), v.size(), v[0].precision());
    for(uint i=0; i<v.size(); i++){
        res[i][i] = v[i];
    }
    return res;
}

vector<Matriz> Matriz::autovalores(TFloat tol, int maxIte){
    vector<TFloat> dant=this->diagonal();
    vector<TFloat> d(alto);
    Matriz Qac(alto);
    Matriz A1 (alto, ancho, precision);
    A1.matriz = this->matriz;
    vector<Matriz> QR (2);
    for (int i = 0; /*((norma1(Qac.restar(dant, d))).dbl() > tol) && */(i < maxIte); i++){
        cout << "Voy por la iteración " << i << " del cálculo de autovectores" << endl;
        QR = A1.qr_rapido();
        A1 = QR[1].multiplicar(QR[0]);
        Qac = Qac.multiplicar(QR[0]);
        dant = d;
        d = A1.diagonal();
        cout<<i<<endl;
    }
    Matriz dRes(alto, ancho, precision);
    Matriz vRes(alto, ancho, precision);
    dRes = dRes.diagonal(d);
    vector<Matriz> res;
    res.push_back(dRes);
    res.push_back(Qac);
    res[0].mostrar(cout);
    res[1].mostrar(cout);
    return res;
}

void Matriz::operator= (const Matriz& otra){
    alto=otra.alto;
	ancho=otra.ancho;
    precision=otra.precision;
    matriz=otra.matriz;
}

vector< vector<TFloat> > Matriz::dameColumnas(){
    vector< vector<TFloat> > columnas;
    vector<TFloat> tmp;
    for(int j=0; j<ancho; j++){
        for (int i = 0; i < alto; i++){
        tmp.push_back(matriz[i][j]);
        }
        columnas.push_back(tmp);
        tmp.clear();
    }
    return columnas;
}

TFloat norma1(vector<TFloat> v){
    float acum = 0;
    for (uint i=0; i<v.size(); i++){
        acum = acum + v[i].dbl();
    }
    TFloat t(acum, v[0].precision());
    return t;
}

TFloat norma2(vector<TFloat> v){
    int tam = v.size();
    float acum = 0;
    for (int i=0; i<tam; i++){
        acum = acum + v[i].dbl()*v[i].dbl();
    }
    TFloat t(sqrt(acum), v[0].precision());
    return t;
}

vector<Matriz> Matriz::qr(){
    vector<TFloat> c;
    TFloat alpha(precision);
    TFloat converter(precision);
    Matriz Aux(alto, ancho, precision);
    Matriz Q(alto);
    vector<Matriz> res;

    Matriz R(alto, ancho, precision);
    R.matriz = this -> matriz;

    //float converter2;

    for(int i=0; i <alto-1 ; i++) {
        cout << "entré al for de qr" << endl;

        c = R.dameColumnas()[i];

        converter = 0;

        for (int j=0; j<=i-1; j++)
            c[j]=converter;

        alpha =norma2(c);

        converter = -1;

        if (c[i].dbl() >=0)
            alpha = converter*alpha;

        c[i] = c[i]+alpha;

        converter = norma2(c);

        for (uint j=0; j<c.size(); j++)
            c[j]=c[j]/converter;

        cout << "terminé el segundo for dentro de qr" << endl;

        Matriz Id(alto);
        Matriz v(c);

        converter = 2;

        Aux = v.multiplicar(v.trasponer()).multiplicar(converter);

        cout << "primera multiplicación" << endl;

        Id = Id.restar(Aux);


        cout << "primera resta" << endl;

        Matriz Idt = Id.trasponer();

        cout << "trasposición" << endl;

        Q = Q.multiplicar(Idt);

        cout << "segunda multiplicación" << endl;

        R = Id.multiplicar(R);

        cout << "tercera multiplicación" << endl;

        cout << "terminé las multiplicaciones" << endl;

        //Q = I-2*v*v'

        //Matriz show(c);
        //show.mostrar(cout);
    }

    res.push_back(Q);
    res.push_back(R);
/*
    //Q.mostrar(cout);
    //R.mostrar(cout);
*/
    return res;
}

int signo(TFloat t){
    if (t.dbl() < 0) return -1;
    else return 1;
}

void normalizar(vector<TFloat> &v){
    TFloat norm = norma2(v);
    for (uint i = 0; i < v.size(); i++)
        v[i] = v[i] / norm;
}

vector<Matriz> Matriz::qr_rapido(){
    Matriz Q(alto);     //inicializo Q como la identidad
    Matriz R(*this);    //inicializo R como copia de mí mismo
    vector<Matriz> Qs;  //vector para calcular las Qs al final

    for (int k = 0; k < alto - 1; k++){
        cout << "Voy por la iteración " << k << " del cálculo de la factorización QR" << endl;
        //defino el vector de R donde quiero poner ceros
        vector<TFloat> v(alto - k);
        //copio la columna de mi matriz
        for (uint i = 0; i < v.size(); i++)
            v[i] = R[i+k][k];
        //modifico la primera componente
        v[0] = v[0] + signo(v[0])*norma2(v).dbl();
        //normalizo v
        normalizar(v);
        //actualizo la matriz R
        Matriz Aux(R);   //para hacer las cuentas modificando R al mismo tiempo
        for (int i = k; i < R.alto; i++)
            for (int j = k; j < R.ancho; j++){
                //calculo (v.v'.R)[i][j]
                TFloat acum(0, precision);
                for (uint z = 0; z < v.size(); z++){
                    acum = acum + v[i-k]*v[z]*Aux[z+k][j];
                }
                //Le resto 2*(v.v'.R)[i][j] a R[i][j]
                R[i][j] = R[i][j] - 2*acum.dbl();
            }
        //R.mostrar(cout);
        //guardo la Q_k que usé
        //Qs.push_back(Q);
    }
    //TO-DO: Calculo Q

    vector<Matriz> res;
    res.push_back(Q);
    res.push_back(R);

    //Q.mostrar(cout);
    R.mostrar(cout);

    return res;
}
